import ps1
import re
import sys
import itertools

Seq = ps1.readSeqsFromFile(sys.argv[1]);


    
p = 0.98; q = 0.999; nSymb = 8;
    
T = [[0.180*p, 0.274*p, 0.426*p, 0.120*p, (1-p)/4, (1-p)/4, (1-p)/4, (1-p)/4],
     [0.171*p, 0.368*p, 0.274*p, 0.188*p, (1-p)/4, (1-p)/4, (1-p)/4, (1-p)/4],
     [0.161*p, 0.339*p, 0.375*p, 0.125*p, (1-p)/4, (1-p)/4, (1-p)/4, (1-p)/4],
     [0.079*p, 0.355*p, 0.384*p, 0.182*p, (1-p)/4, (1-p)/4, (1-p)/4, (1-p)/4],
     [(1-q)/4, (1-q)/4, (1-q)/4, (1-q)/4, 0.300*q, 0.205*q, 0.285*q, 0.210*q],
     [(1-q)/4, (1-q)/4, (1-q)/4, (1-q)/4, 0.322*q, 0.298*q, 0.078*q, 0.302*q],
     [(1-q)/4, (1-q)/4, (1-q)/4, (1-q)/4, 0.248*q, 0.246*q, 0.298*q, 0.208*q],
     [(1-q)/4, (1-q)/4, (1-q)/4, (1-q)/4, 0.177*q, 0.239*q, 0.292*q, 0.292*q]]
     
def argmax(array):
    max = array[0]; index = 0;
    for i in range(1,len(array)):
        if(array[i]>max): 
            max = array[i]; index = i;
        
    return index

# akl: A+ C+ G+ T+ A- C- G- T
def E(i,c):
    
    if(c=='A'):
        if(i==0 or i==4): return 1;
        else: return 0;
        
    elif(c=='C'):
        if(i==1 or i==5): return 1;
        else: return 0;
    
    elif(c=='G'):
        if(i==2 or i==6): return 1;
        else: return 0;
        
    elif(c=='T'):
        if(i==3 or i==7): return 1;
        else: return 0;

    
#E = [[1, 0, 0, 0, 1, 0, 0, 0],
#     [0, 1, 0, 0, 0, 1, 0, 0],
#     [0, 0, 1, 0, 0, 0, 1, 0],
#     [0, 0, 0, 1, 0, 0, 0, 1]]
# Viterbi
    V = list(); StateViterbi = list();
    V = [[0,0,0,0,0,0,0,0]]; 
    for k in range(0,nSymb):
        V[0][k] = E(k,Seq[iseq][0])*1/8;
    
    StateViterbi.append(argmax(V[0]));

    for i in range(0,len(Seq[iseq])-1):
        V.append([0,0,0,0,0,0,0,0]);
        temp = [0,0,0,0,0,0,0,0];
        for l in range(0,nSymb):
            for k in range(0,nSymb):
                temp[k] = V[i][k]*T[k][l];
        #for k in range(0,nSymb): temp[k] = temp[k]/sumtemp
    
            V[i+1][l] = E(l,Seq[iseq][i+1])*max(temp);
        sumV = sum(V[i+1]);
        for l in range(0,nSymb):
            V[i+1][l] = V[i+1][l]/sumV    
        
        StateViterbi.append(argmax(V[i+1]));

    del V; del StateViterbi;


for iseq in range(0,len(Seq)):
    
#Viterbi    
    V = list(); 
    V = [[0,0,0,0,0,0,0,0]]; 
    for k in range(0,nSymb):
        V[0][k] = E(k,Seq[iseq][0])*1/8;
    

    for i in range(0,len(Seq[iseq])-1):
        V.append([0,0,0,0,0,0,0,0]);
        temp = [0,0,0,0,0,0,0,0];
        for l in range(0,nSymb):
            for k in range(0,nSymb):
                temp[k] = V[i][k]*T[k][l];
        #for k in range(0,nSymb): temp[k] = temp[k]/sumtemp
    
            V[i+1][l] = E(l,Seq[iseq][i+1])*max(temp);
        sumV = sum(V[i+1]);
        for l in range(0,nSymb):
            V[i+1][l] = V[i+1][l]/sumV    
        

    
    
    Z = list(); 
    for i in range(len(Seq[iseq])-1,-1,-1): Z.append(-1);
    
    Z[len(Seq[iseq])-1] = argmax(V[len(Seq[iseq])-1]);
    # BACKWARD
    for i in range(len(Seq[iseq])-2,-1,-1):
        for k in range(0,nSymb):
            temp[k] = E(Z[i+1],Seq[iseq][i+1])*V[i][k]*T[k][Z[i+1]];
        Z[i]=argmax(temp);
        
    del V;
        
# FORWARD PASS
    Alpha = list(); Alpha = [[0,0,0,0,0,0,0,0]];
    for k in range(0,nSymb):
        Alpha[0][k] = E(k,Seq[iseq][0])*1/8;

    for i in range(0,len(Seq[iseq])-1):
        Alpha.append([0,0,0,0,0,0,0,0]);
        temp = [0,0,0,0,0,0,0,0];
        for l in range(0,nSymb):
            for k in range(0,nSymb):
                temp[k] = Alpha[i][k]*T[k][l];
        #for k in range(0,nSymb): temp[k] = temp[k]/sumtemp
    
            Alpha[i+1][l] = E(l,Seq[iseq][i+1])*sum(temp);
            sumF = sum(Alpha[i+1]);
        for l in range(0,nSymb):
            Alpha[i+1][l] = Alpha[i+1][l]/sumF    


        
# BACKWARD PASS
    Beta = list(); 
    for i in range(0,len(Seq[iseq])):
        Beta.append([1,1,1,1,1,1,1,1,]);
    

    for i in range(len(Seq[iseq])-2,-1,-1):
        for k in range(0,nSymb):
            temp = [0,0,0,0,0,0,0,0];
            for l in range(0,nSymb):
                temp[l] = T[k][l]*E(l,Seq[iseq][i+1])*Beta[i+1][l];
            Beta[i][k] = sum(temp);
    
        sumB = sum(Beta[i]);
        for k in range(0,nSymb):
            Beta[i][k] = Beta[i][k]/sumB;
            
    P = list();   
    for i in range(0,len(Seq[iseq])):
        P.append([1,1,1,1,1,1,1,1,]);    
        for k in range(0,nSymb):
            P[i][k] = Alpha[i][k]*Beta[i][k];
        sumP = sum(P[i]);
        for k in range(0,nSymb):
            P[i][k] = P[i][k]/sumP;
     
    del Alpha; del Beta;
    State = list(); Differ = list();
    for i in range(0,len(Seq[iseq])):
        State.append(argmax(P[i])); 
    for i in range(0,len(Seq[iseq])):
        if(State[i] !=Z[i]):
            Differ.append(i); 
    
        
    print(); print(); print(); print(); print(); print(); print("Gene :" + str(iseq)); 
    print(Seq[iseq]);
    print(); print("Viterbi Path: " + str(Z));
    print(); print("Posterior Path: " + str(State));
    print(); print("Positions where Viterbi & Posterior Decoding Differ");
    print(Differ);
    print("Total Difference: " + str(len(Differ)));
    
    
        
    del P; del State;  del Z; del Differ;